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AN  IMPROVED  TOEPLITZ  APPROXIMATION  METHOD1 


K.  S.  Arun 

Electrical  and  Computer  Engineering 
Coordinated  Science  Laboratory 
Univ.  of  Illinois  at  Urbana-Champaign 
Urbana,  IL-61801 

ABSTRACT 

In  this  paper,  we  suggest  a  modification  of  the  Toe- 
plltz  approximation  method  for  estimating  frequencies  of 
multiple  sinusoids  from  covariance  measurements.  The 
method  constructs  a  state-feedback  matrix  following  a 
low-rank  approximation  of  the  Toeplitz  covariance  matrix 
via  singular  value  decomposition.  Ideally,  the  eigenvalues 
of  this  state-feedback  matrix  will  be  on  the  unit  circle  in 
the  complex  plane,  and  the  angles  that  they  make  with  the 
real  axis  will  be  equal  to  the  unknown  sinusoid  frequen¬ 
cies.  The  modification  proposed  here  exploits  this  prior 
knowledge  of  the  modulus  of  the  eigenvalues,  and  guaran¬ 
tees  that  even  in  the  presence  of  noise,  the  eigenvalues  of 
the  estimated  state-feedback  matrix  will  lie  on  the  unit 
circle. 


1.  INTRODUCTION 

The  problem  of  retrieving  multiple  sinusoids  (with 
frequencies  close  to  each  other)  from  perturbed  time-series 
or  covariance  information  is  of  special  interest  in  a  vast 
range  of  signal-processing  applications.  Very  often  the 
covariance  sequence  may  have  to  be  estimated  from  time- 
series  data,  as  in  Doppler  processing  in  radar.  It  is  nol 
uncommon,  however,  to  encounter  applications  in  which 
the  (time-series)  data  are  not  measurable  while  the  covari¬ 
ance  information  is  directly  available.  Such  situations 
arise  in  astronomical  star  bearing  estimation,  interference 
spectroscopy,  and  some  sensor  array  applications. 

In  recenl  years,  there  has  been  a  great  deal  of  interest 
in  model-based  sinusoid  retrieval.  Models  convert  the 
non-linear  problem  of  estimating  Ihe  sinusoid  frequencies 
into  a  simpler  problem  of  estimating  the  parameters  of  a 
linear  model  (l).  The  second  step  in  all  model-based 
methods  is  the  extraction  of  the  desired  information  (the 
frequencies)  from  the  estimated  model  parameters  (2). 

Both  steps  are  important  for  the  overall  success  of  a 
model-based  melhod.  111-conditioning  at  either  step  can 
adversely  effect  the  overall  performance  of  the  method 
and  should  be  avoided.  The  reliability  of  the  first  step 
depends  on  the  estimation  procedure,  and  that  of  the 
second  step  on  the  sensitivity  of  the  desir"d  information  to 
the  model  parameters  [3].  A  popular  model  for  the 
sinusoid  retrieval  problem  is  the  linear  prediction  model 
6rst  used  by  Prony  in  1 88 1 . 
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y(t)=  £  aky(t-k) 

k=l 

whose  parameters  may  be  reliably  estimated  by  the 
method  of  Tufts  and  Kumaresan  [4],  The  roots  of  the 
polynomial  formed  from  these  parameters  are  ideally 
expected  to  be  on  the  unit  circle  in  the  complex  plane,  and 
the  angles  that  they  make  with  the  real  axis  should  equal 
the  sinusoid  frequencies. 

2.  STATE-SPACE  REPRESENTATION 
It  turns  out  that  the  sinusoidal  model  is  a  very  spe¬ 
cial  case  of  the  general  linear  rational  model,  and  that  just 
as  there  are  alternate  parameterizations  of  linear  systems, 
there  also  are  alternate  parameter  sets  Tor  the  sinusoidal 
model  as  well.  Just  as  there  is  a  state-space  representation 
for  every  realization  of  a  linear,  rational  system,  there  is 
also  a  state-space  representation  for  every  realization  of 
the  sinusoidal  model.  The  stale-space  representation  of  the 
special  model  for  sinusoidal  signals  (frequencies:  o>,. 
l-  1 .2 . n )  is: 

x(k  +  1 )  =  Fx(k) 
y(k)  =  hx(k) ' 

where  the  order  of  the  model  p  is  twice  the  number  ol 
sinusoids,  and  the  eigenvalues  of  F  are  of  unit  magnitude 

and  equal  c1,l‘\  1-1.2 . n.  The  sinusoidal  signal  y(t)  is  the 

model's  zero-input  response  to  some  non-zero  initial  condi¬ 
tion  x(0).  In  fact,  we  have 

y(t)  =  hF’x(O).  t>0. 

and  the  covariance  r(m)  of  the  sinusoidal  signal  satisfies 

r(m)  =  hFmPh'  m>0  (1) 

where  P  is  the  state-covariance  matrix,  and  the  superscript 
t  denotes  the  Hermilian  transpose. 

The  linear  prediction  model  is  a  canonical  realization 
of  the  above,  with 

x(t)  =  j  y( t—  1 )  y(t— 2)  ■  •  •  y(t-p)  j'. 

a,  a2 
1  0 
0  1 


0  0 
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Another  canonical  realisation  which  avoids  the  second  step 
completely  is 

F  =  diag  (e2'1'1).  h  =  (1,1 . l). 

However,  the  first  step  of  estimating  the  model  parameters 
becomes  difficult  for  this  realization.  This  non-uniqueness 
of  the  parameter  triple  (F.x(0),h)  that  characterizes  the 
realization,  allows  one  to  choose  a  realization  that  makes 
both  steps  of  the  model-based  method  reliable.  The  state 
space  parameters  can  be  estimated  from  covariance  data  by 
using  the  factorizations  derived  below. 

Using  Eq.  (1)  for  the  covariance  lags,  and  noting  that 
the  state  covariance  matrix  P  satisfies  P  =  FPF1.  it  can  be 
shown  that 

r(— m)  =  hF'mPh’.  (2) 

Using  Eqs.  ( 1 )  and  (2).  the  Toeplitz  covariance  matrix 


r(0) 

r(-l) 

r(— 2)  . 

.  r(— n) 

r(l) 

r(0) 

r(— 1 )  . 

.  r(— n+1) 

r(2) 

Kl) 

r(0) 

.  r(— n+2) 

r(n) 

r(n— 1 ) 

r(n— 2)  . 

.  r(0) 

can  be  factorized  as  shown  below. 


h 

hF 

hF2 


jph1  F"'Ph'  F_2Ph’ 


.  .  .  y- 


r  I 

k  =  er. 

Both  0  and  T  are  full  rank,  and  so  the  rank  of  R  is  equal 
to  the  model  order  p,  which  is  twice  the  number  of 
sinusoids  [5).  Observe  that  the  ilh  row  of  ©  is  hF1-1.  and 
the  i’h  column  of  T  is  F_l+,Ph\  so  thai  F  may  be  obtained 
by  solving  the  overdetermined  system  of  equations 

e,F  =  62  (3) 

where  0,  (02)  is  obtained  from  0  by  deleting  the  last 
(first)  row,  or  by  solving 

F-'r,  =r2.  (4) 

where  and  T2  are  defined  in  a  manner  similar  to  0,  and 
02  respectively. 


3.  ORIGINAL  TAM 

The  above  discussion  indicates  how  the  sinusoid  fre¬ 
quencies  may  be  obtained  from  covariance  data,  when  the 
information  is  exact.  Let  the  singular  value  decomposition 
(SVD)  of  the  Toeplitz  matrix  R  (which  is  also  the  eigen- 
decomposition  when  the  covariances  are  exact)  be 
R=UEV\  where  E  contains  only  the  non-zero  singular 
values.  Hence,  the  dimensions  of  diagonal  matrix  E  will 
equal  the  rank  of  R,  which  is  also  the  order  of  the  model  p. 
Different  factorizations 

©  =  UE'^Q,  T  =  Q-'E'^V' 

of  the  Toeplitz  matrix  arc  possible,  and  each  will  lead  to  a 


diUercnl  state-space  realization  of  the  order-p  sinusoidal 
model  as 

f  =  ©,l©2=  (r2r  fr*.  (5) 

Here,  the  superscript  L  (similarly  R)  denotes  any  left- 
inverse  (or  right-inverse).  They  exist  because  the  matrices 
under  consideration  have  full  rank.  The  frequencies  of  the 
sinusoids  may  then  be  found  as  the  angles  of  the  eigen¬ 
values  of  the  F  matrix. 

Theoretically,  the  algorithm  should  work  for  any 
choice  of  Q-  For  instance,  the  choice  Q  =  E'^Vf’,  where  Vf 
denotes  the  matrix  formed  from  the  first  p  rows  of  V,  will 
lead  to  the  canonical  realization  with  the  linear  prediction 
parameters.  The  so-called  balanced  choice  0  =  I  is  the  pre¬ 
ferred  choice  for  reducing  the  sensitivity  of  the  sinusoid 
frequencies  lo  the  model  parameters. 

It  is  well  known  in  the  numerical  analysis  literature 
that  the  eigenvalues  of  a  normal  matrix  are  least  sensitive 
to  perturbations  in  the  matrix  entries  [fa].  Recall  that  F  is 
a  no'rma)  matrix  if  it  satisfies  FF’  =  F’F.  Since  the  state- 
feedback  matrix  of  the  sinusoidal  model  has  eigenvalues  of 
unit  modulus,  this  amounts  to  requiring  that  the  matrix  be 
unitary.  Note  that  for  every  realization  of  the  sinusoidal 
model,  F  has  eigenvalues  on  the  unit  circle,  which  is  a 
necessary  condition  for  it  to  be  unitary.  It  is  not  a 
sufficient  condition  however,  and  not  all  realizations  of  the 
sinusoidal  model  have  unitary  state-feedback  matrices; 
but  it  turns  out  that 


Theorem  I:  The  F  matrix  obtained  from  any  symmetric 

factorization  (0  =  I"’)  of  the  square  covariance  matrix  R  is 
unitary. 

Proof:  Let  0  =  T'.  Then  0,  =  T,'  and  02  =  r2.  Com¬ 

bining  this  with  Eq.  (5),  we  have  F  =  (F-1)'  which  implies 
F  is  unitary.  Q 

The  Toeplitz  approximation  method  (TAM)  of  [7,5] 
exploits  these  facts  to  reliably  estimate  the  sinusoid  fre¬ 
quencies  from  inexact  covariances.  Inexactness  can  be 
caused  in  practice,  by  a  number  of  factors,  additive  noise 
in  the  data,  errors  in  estimating  the  covariances,  finite  pre¬ 
cision  errors,  and  others.  At  first,  TAM  performs  an  SVD 
of  R.  and  retains  the  p  principal  components  (i.e.,  the  p 
largest  singular  values  and  the  corresponding  singular  vec¬ 
tors).  SVD  is  preferred  to  eigendecomposition  because  in 
the  presence  of  perturbations,  R  may  not  be  non-negative 
definite.  Let  the  singular  vectors  and  singular  values  after 
the  low-rank  approximation  be  U.E.  and  V.  Next,  TAM 
picks  0  =  UEV*.  and  looks  for  an  approximate  solution  to 
Eq.  (3).  The  approximation  criterion  used  is  least-squares, 
and  the  TAM  estimate  is 

Ft  am  =  ®\&i-  «>) 

where  the  superscript  t  denotes  the  pseudo-inverse.  The 
sinusoid  frequency  estimates  are  then,  the  angles  of  the 
eigenvalues  of  F,  AM. 


Theorem  1  indicates  that  in  the  noiseless  situation. 

1  AM's  choice  of  factor  0  as  UEV'  will  make  the  F-matrix 
unitary.  In  the  low  SNR  case  as  well,  if  there  is  an  exact 
solution  to  Lq.  (3),  (which  is  the  case  when  no  low-rank 
approximation  is  made,  R  is  non-negative  definite  and 
exactly  equal  lo  UEV1,  lor  instance)  then  the  solution  will 

Cedes 
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be  unitary.  But  whenever  a  principal-components  approx¬ 
imation  is  made  and  there  is  no  exact  solution  to  l:q.  (3). 
the  least-squares  solution  FTAV1  will  gcnerically  not  be 
unitary. 

4.  IMPROVED  TAM 

The  modification  proposed  in  this  section  is  the  inclu¬ 
sion  of  the  unitary  constraint  on  F  in  the  least-squares 
solution  of  Eq.  (3).  Such  a  constrained  least -squares 
problem  also  arises  in  motion  estimation  and  a  simple 
analytical  solution  was  derived  in  (8).  We  include  the 
problem  and  the  solution  here,  but  refer  the  reader  to  [8] 
for  the  derivation  of  the  algorithm. 

Problem:  Given  two  real-valued  n  by  p  matrices  ©!  and 

©2,  n^p.  find  p  by  p  matrix  F  that  minimizes 

Trace  j(©,F  -  02)’(©,F  -  02)  J. 

subject  to  the  constraint  that  F  be  unitary. 

Solution.  Calculate  the  p  by  p  inauix 

H  =  02’e,. 

Let  the  SVD  of  H  be  UHIHV^,  where  IH  includes  all  p 
singular  values,  even  if  they  are  zero.  Then  the  solution  to 
the  constrained  least-squares  problem  is 

Funiury  =  V„U,5.  (?) 

We  propose  the  use  of  this  solution  in  place  of  Eq.  (6) 
to  guarantee  that  the  estimated  F-matrix  is  unitary  and 
consequently  does  have  eigenvalues  on  the  unit  circle.  To 
reduce  the  effects  of  noise  on  the  estimale  of  F,  TAM's 
choice  of  0  is  also  slightly  modified  as  follows.  The  com¬ 
plete  algorithm  is 

Improved  TAM  algorithm: 

1.  Compute  the  SVD  of  R 

r 

R  =  E  un°Kvk' 

k=l 

where  the  rank  r  is  at  least  as  large  as  p.  and  the  singular 
values  are  numbered  in  order  of  decreasing  magnitude. 
Then,  denote 

U=  [u,  -  •  ■  upj.  I  =  diag(ok  -  op+1;  k=l . p); 

and  choose 

0=  UI'" 

2.  Let  the  SVD  of  H  =  ©,'©,  be  U„E„V,;.  where  E„ 
includes  all  p  singular  values.  Then  the  new  estimale  ol 
the  state-feedback  matrix  is 

^  unitary  =  VI|U||- 

A  heuristic  reasoning  for  subtracting  the  (p+l)-ih 
singular  value  in  the  construction  of  0  is  the  following. 
When  there  is  additive  white  noise  in  the  data  and  the 
covariances  arc  exact,  all  singular  values  ol  K  gel 
translated  up  by  the  noise  variance,  but  the  singular  vec¬ 
tors  are  not  effected.  So.  if  0  Is  constructed  from  the  p 
principal  components,  then,  Eq.  (3)  will  have  an  exact 
solution.  But  if  op+j  is  not  subtracted  off,  the  solution 
will  not  be  unitary. 


5.  SIMULATIONS 

Our  simulations  have  shown  noticeable  improvement 
in  TAM's  performance  when  exact  covariances  are  avail¬ 
able;  but  negligible  improvement  when  covariances  have  to 
be  estimated  from  low  SNR  data. 

Example  I.  In  the  first  example,  it  is  assumed  that 
exact  covariances  are  available,  corresponding  to  a  single 
sinusoid  in  an  additive  AR(  1 )  process. 

r(m)  =  cos(2?rfm)  +  p(0.5)lml 

The  improved  TAM  and  the  original  TAM  were  used  on 

covariances  r(0) .  r(I2);  corresponding  lo  different 

values  of  the  sinusoid  frequency  f.  The  Toeplitz  matrix 
used  was  of  size  13  x  13,  and  a  rank-2  approximation  was 
employed  in  each  algorithm.  The  simulations  were  per¬ 
formed  in  double  precision  on  a  VAX  780  using  FORTRAN 
77.  The  1MSL  routines  LSVDF,  and  E1GRF  were  used  in 
the  programs. 

Table  I  gives  the  resulis  for  SNR  10  dB  (p  =  0.1). 
Our  results  indicate  that  TAM  fails  to  detect  the  sinusoid 
at  frequencies  around  10--’  (i.e.,  both  eigenvalues  were 
real),  and  the  improved  TAM  is  able  to  resolve  the  two 
complex  exponentials  at  much  smaller  frequencies,  bui 
appears  to  make  large  errors  in  the  frequency  estimate. 

Table  II  gives  the  results  for  low  SNR  for  a  fixed  fre¬ 
quency  f-0.2.  Both  methods  appear  to  fail  ai  the  same 
SNR  in  this  simulation. 

Example  2.  In  this  example,  data  corresponding  to  two 
cqui-amplitude  real  sinusoids  in  additive  white  noise  was 
used  as  the  starling  point  for  the  estimation  schemes. 

y(n)=  Acos(2nT,n  )  +  Acos(2rrf2n+</>)  +  w(n), 

where  <b  is  a  random  variable  uniformly  distributed  over 
[— tr.+7rj,  and  w(n)  is  /cro-mcan  Gaussian  white  noise, 
independent  of  <t>-  Only  48  consecutive  observations  from 
a  single  sample  sequence  arc  given;  and  the  covariance  lags 
have  to  be  estimated  by  temporal  averaging  over  the  single 
record.  For  the  experiment.  IMSI  routines  GGUBS  and 
GGNML  were  used  to  generate  <b  and  w(n)  respectively. 
The  original  TAM  as  well  as  the  improved  TAM  were  both 
employed  to  estimate  the  sinusoid  frequencies  from 
unbiased  covariance  estimates,  using  rank-2  approxima¬ 
tions.  The  experiment  was  repeated  50  times  using 
independent  realizations  of  <t>  and  w(n).  and  ihe  arithmetic 
mean  and  standard  deviation  was  computed  for  the  50 
estimates  of  each  frequency.  The  simulations  were  per¬ 
formed  in  single  precision  on  a  60-bil  CYBER  175  using 
FORTRAN  77. 

Table  III  gives  the  results  for  SNR  10  dB.  f,  =0.125. 
f2  =  0.135.  and  matrix  size  32  by  32.  Fable  IV  gives  the 
results  for  SNR  3  dB.  ft  =0.125,  f2  =  0.145,  and  matrix 
size  24  by  24.  The  results  seem  to  indicate  that  any 
improvement  over  ihe  original  TAM  is  marginal  in  this 
example.  A  possible  explanation  is  that  TAM  has  achieved 
(or  is  very  close  to)  the  CTamcr-Ruo  lower  bound  and  can¬ 
not  be  further  improved.  Another  possibility  is  that  the 
state-feedback  matrix  estimated  by  the  the  original  TAM 
is  already  close  to  unitary,  even  before  imposing  the  con¬ 
straint.  Both  possibilities  arc  the  subject  of  current  inves¬ 
tigation. 
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True  f 

0.20000 

0.02000 

0.00200 

0.00100 

0.00050 

0.00020 

E2J5HSZEBI 

0.19979 

0.02004 

0.00184 

0.00000 

0.00000 

0.00000 

0.19980 

0.01988 

0.00206 

0.00158 

0.00146 

0.00143 

Table  1  (SNR  -  10  dB) 


SNR 

10  dB 

0  dB 

-2  dB 

-4  dB 

-5  dB 

<EHn 

0.19979 

0.19767 

0.19583 

0.18968 

0.16661 

Improved  TAM  estimate 

0.19980 

0.19781 

0.19607 

0.19000 

0.16362 

Table  II  (f  -0.2) 


- - 

r? 

Mean 

0.12848 

0.13761 

f OTtfatfSi 

Std.  dev. 

m-rnrei 

5  35  x  10-6 

Mean 

0.12854 

0.13743 

Std.  dev. 

tftinra 

5.21  x  lO-6 

Mean 

0.12713 

0.16271 

Std.  dev. 

4.34  x  10-5 

7.42  x  10-3 

Mean 

0.12439 

0.15275 

F.HWliTVli 

Std.  dev. 

3.64  x  10-" 

3.71  x  10~3 

Table  111  (f,=0.125.  f2-0.135)  Table  ,v  (f1=o.i25.  f2=0.145) 
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6.  SUMMARY 

In  summary,  this  paper  examined  the  imposition  of  a 
unitary  constraint  on  the  state-feedback  matrix  con¬ 
structed  by  the  Toeplitz  approximation  method  for 
estimating  ihe  frequencies  of  multiple  sinusoids.  The  con¬ 
straint  is  a  means  of  incorporating  given  knowledge  about 
the  sinusoidal  model,  namely  that  the  eigenvalues  of  the 
state-feedback  matrix  are  on  the  unit  circle.  This  informa¬ 
tion  can  also  be  partially  incorporated  in  linear-prediction 
methods  for  the  same  problem.  It  has  been  previously 
suggested  that  a  symmetry  constraint  be  imposed  on  the 
estimated  linear  prediction  polynomial,  because  it  is 
known  that  the  polynomial  must  necessarily  be  symmetric 
to  have  roots  on  the  unit  circle.  However,  this  symmetry 
condition  is  only  necessary  and  not  sufficient  for  ensuring 
that  the  roots  have  unit  modulus.  For  instance,  a  sym¬ 
metric  polynomial  may  have  roots  at  <*  and  -L  instead  or 
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unit-modulus  roots.  In  contrast,  the  constraint  studied 
here  is  both  necessary  and  suBicient  as  demonstrated  by 
Theorem  1, 

Preliminary  numerical  simulations  seem  to  indicate 
that  an  improvement  is  possible,  when  covariances  are 
directly  measured.  In  every  example,  the  modified  TAM 
was  at  least  as  good  as  the  original  TAM,  both  in  resolu¬ 
tion  and  in  estimation  error.  The  modification  proposed 
here  may  be  particularly  useful  in  applications  where  the 
covariances  are  obtained  by  (pseudo-)  ensemble  averaging 
instead  of  temporal  averaging.  Further  experimentation 
and  analysis  Is  needed  to  determine  the  kinds  of  problems 
for  which  the  extra  computational  expense  pays  off  in 
marked  improvement. 
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